Regressão Linear com Python¶
Autor: Ronisson Lucas Calmon da Conceição
1. Estimadores de MQO¶
Características dos estimadores de MQO:
- Fácil aplicação
- Simples interpretação e propriedades desejáveis
1.1 Regressão Linear Simples¶
Introdução
- $\mathbf{y}$: representa o fenômeno/variável target que queremos analisar em função de uma variável explicativa $\mathbf{x}$ (ou mesmo conjunto de variáveis explicativas) e um termo de perturbação $\mathbf{u}$ (ou termo de erro): $$ y = f(x)+u $$
- Vamos focar inicialmente no problema que busca analisar a relação entre duas variáveis descrita por um modelo de regressão linear simples:
$$ y = \beta_0+\beta_1 x +u$$
Note que: $\mathbf{y}$ é a variável dependente e $\mathbf{x}$ é a variável explicativa (que prova as variações em $y$).
Exemplos - $\text{salário(w)} = \beta_0+\beta_1\text{educação}+u$: o parâmetro $\beta_1$ refere-se ao efeito do grau de escolaridade do indivíduo no salário.
- $\text{taxa_crime(e)} = \beta_0+\beta_1\text{quantidade_polícia}+u$: o parâmetro $\beta_1$ refere-se a redução da taxa de criminalidade em função de um aumento na quantidade de força polícial.
- $\text{vendas(s)} = \beta_0+\beta_1\text{temperatura}+u$: para algum tipo de negócio o aumento da temperatura pode ter um efeito positivo na quantidade de vendas.
Terminologia do modelo de regressão linear simples:
| y | x |
|---|---|
| Variável Dependente | Variável Independente |
| Variável Explicada | Variável Explicativa |
| Variável de Resposta | Variáve de Controle |
| Variável Prevista | Variável Previsora |
| Regressando | Regressor |
Ou ainda:
$\underbrace{y}_{\text{target}} = \underbrace{\beta_0+\beta_1}_{\text{coeficientes}} \underbrace{ x}_{\text{ input}} +\underbrace{u}_{\text{termo de erro}}$
O termo de erro $\mathbf{u}$ diz respeito a outros fatores que afetam $\mathbf{y}$ mas não estão explícitos na equação da regressão.
Relação linear entre $\mathbf{y}$ e $\mathbf{x}$: $\mathbf{x}$ altera a variável $\mathbf{y}$ linearmente. Uma implicação desta premissa seria que independentemente do valor inicial de $\mathbf{x}$ teremos sempre o mesmo impacto de desta variável sobre $\mathbf{y}$.
Lembre que a relação é linear nos parâmetros do modelo.
Considere um modelo populacional $y = \beta_0+\beta_1 x+u$ e uma amostra desta população $(x_i,y_i), i = 1, \dots , n$. Podemos obter que os estimadores de MQO para uma regressão linear simples:
$\min_{\hat{\beta}_0, \hat{\beta}_1} \sum_{i=1}^{n}\hat{u}_i^2 = \sum_{i=1}^{n}\left(y_i-\hat{\beta}_0-\hat{\beta}_1x_i\right)^2$
$\dfrac{\partial \sum_{i=1}^{n}\hat{u}_i^2}{\partial \hat{\beta}_0} = 0$
$\dfrac{\partial \sum_{i=1}^{n}\hat{u}_i^2}{\partial \hat{\beta}_0} = 0$
Modelo estimado:
Resíduo do modelo:
Assim o resíduo da regressão é definido como a diferença entre o valor $y_i$ (valor real) e o valor previsto/ajustado $\hat{y}_i$. Se $\hat{u}_i>0 \Rightarrow$ então a reta subestima $y_i$; por outro lado, se $\hat{u}_i<0$ então a reta superestima $y_i$. O mundo ideal seria $\hat{u}_i = 0$, mas na prática teremos para a maior parte dos casos $\hat{u}_i\neq 0$.
Hands On!¶
Nosso objetivo é estimar um modelo para explicar o salário com base no tempo de experiência.
import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
import seaborn as sns
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import warnings
warnings.filterwarnings("ignore")
df = pd.read_csv('data.csv')
df.head()
df.describe()
sns.pairplot(df);
# correlação linear de Pearson
df.corr()
figure = go.Figure(data = go.Scatter(x = df['tempo_experiencia'],
y = df['salario'],
mode = 'markers'
),
layout = go.Layout(xaxis = dict(title = 'Tempo de experiência'),
yaxis = dict(title = 'Salário'),
title = 'Relação entre salário e tempo de experiência'
)
)
figure.update_yaxes(tickprefix = '$')
figure
Nosso objetivo é então obter uma reta que melhor se ajuste aos dados.
modelo = smf.ols(formula = 'salario ~ tempo_experiencia', data = df).fit()
modelo.params
# previsão
fitted_values = modelo.fittedvalues
figure = go.Figure(
data = [
go.Scatter(
x = df['tempo_experiencia'],
y = fitted_values,
name = 'Reta ajustada',
mode = 'lines',
line = dict(color = 'black')
),
go.Scatter(
x = df['tempo_experiencia'],
y = df['salario'],
mode = 'markers',
name = 'Amostra'
)
],
layout = go.Layout(
xaxis = dict(title = 'Tempo de experiência'),
yaxis = dict(title = 'Salário'),
title = 'OLS: salário ~ tempo de experiência'
)
)
figure.update_yaxes(tickprefix = '$')
figure
Regressão linear múltipla¶
Na prática teremos mais de uma variável explicativa e podemos então generalizar para um modelo com $K$ variáveis explicativas e incorporar mais mais fatores que podem explicar $y$.
$$y = \beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3+\dots+\beta_Kx_k+u$$Nota: $\beta_0$ é o intercepto do modelo, $\beta_1$ é o parâmetro associado a $x_1$, $\beta_2$ é o parâmetro associado a variável $x_2$, e assim sucessivamente.
Podemos usar notação matricial:
Sendo que:
$ \mathbf {y} = \left[ \begin{array}{c c c} y_{1} \\ y_{2} \\ \vdots\\ y_{n}\\ \end{array}\right] $, $ \mathbf {\mathbf{\beta}} = \left[ \begin{array}{c c c} \beta_{0} \\ \beta_{1} \\ \vdots\\ \beta_{k}\\ \end{array}\right] $, $ \mathbf {u} = \left[ \begin{array}{c c c} u_{0} \\ u_{1} \\ \vdots\\ u_{n}\\ \end{array}\right] $, $ \mathbf {X} = \left[ \begin{array}{c c c} 1 & x_{11}& \ldots& x_{1k}\\ 1 & x_{21}& \ldots& x_{2k}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n1} &\ldots & x_{nk} \end{array}\right] $
Nosso problema de otimização e obtenção dos parâmetros consiste em:
$$\min_{\hat{\beta}}\mathbf{\hat{u'}\hat{u}} = \left(\mathbf{y}-\mathbf{X\hat{\beta}}\right)'\left(\mathbf{y}-\mathbf{X\hat{\beta}}\right)$$Que resolvendo obtemos:
$$\therefore \mathbf {\hat\beta} = (\mathbf {X^{'}}\mathbf {X})^{-1} \mathbf {X^{'}}y $$import statsmodels.api as sm
from sklearn.model_selection import train_test_split
import statsmodels.stats.api as sms
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
def line_plot(xaxis, yaxis, data_frame, title = '',
figsize = (12, 6),
xlabel = '',
ylabel = '',
line_color = 'k'
):
fig, ax = plt.subplots(figsize = figsize)
ax.plot(data_frame[xaxis], data_frame[yaxis],
color = line_color
)
ax.set(title = title, xlabel = xlabel, ylabel = ylabel)
def heatmap_plot(data_frame,
correlation_type = 'pearson',
title = '',
annot = True,
cmap = 'Blues',
figsize = (8, 8)
):
corr = data_frame.corr(method = correlation_type)
fig, ax = plt.subplots(figsize = figsize)
sns.heatmap(corr, ax = ax,
annot = annot,
cmap = cmap
)
ax.set(title = title)
def count_plot(data, labelsize = 12,
xlabel = '',
ylabel = 'Frequência',
style = 'white',
):
sns.set(rc = {
'figure.figsize':(10, 5),
'xtick.labelsize':labelsize,
'ytick.labelsize':labelsize,
})
sns.set_style(style)
count_plot = sns.countplot(x = data)
count_plot.set_title('Countplot de '+data.name, fontsize = 15)
count_plot.set(xlabel = xlabel,
ylabel = ylabel
)
#def distribution_plotter(data, label, bins=None):
def distribution_plotter(data, label, bin_width=200):
sns.set(rc={"figure.figsize": (10, 7)})
sns.set_style("white")
#dist = sns.distplot(data, bins= bins, hist_kws={'alpha':0.2}, kde_kws={'linewidth':5})
dist = sns.histplot(data,
stat = 'density', kde = True,
line_kws = {'linewidth':5},
binwidth = bin_width)
dist.set_title('Distribuição de ' + label + '\n', fontsize=16)
Hands On!¶
df = pd.read_csv('Consumo_cerveja.csv')
df.head()
df.shape
df.drop_duplicates(inplace = True, keep = False)
df.shape
df.isnull().sum()
df.info()
df.dtypes
for column in df.select_dtypes(include = object).columns[1:]:
df[column] = pd.to_numeric(df[column].str.replace(',','.'))
df['Data'] = pd.to_datetime(df['Data'], format = '%Y-%m-%d')
df.columns
df.rename(
columns = {
'Temperatura Media (C)': 'Temperatura_mean',
'Temperatura Maxima (C)': 'Temperatura_Max',
'Temperatura Minima (C)': 'Temperatura_min',
'Precipitacao (mm)': 'Precipitacao',
'Final de Semana': 'Final_semana',
'Consumo de cerveja (litros)': 'Consumo_cerveja'
},
inplace = True
)
df.info()
df.describe()
cols_temperatura = df.iloc[:,1:4]
fig, ax = plt.subplots(figsize = (12, 6))
for column in cols_temperatura:
ax.plot(df['Data'], df[column])
ax.set(title = 'Temperaturas máxima, média e mínima',
xlabel = 'Data'
);
traces = list()
for column in cols_temperatura:
trace = go.Scatter(
x = df['Data'],
y = df[column],
name = column,
mode = 'lines'
)
traces.append(trace)
layout = go.Layout(
xaxis = dict(title = 'Data'),
yaxis = dict(title = 'Temperatura'),
title = 'Temperaturas máxima, média e mínima'
)
fig = go.Figure(data = traces, layout = layout)
fig
line_plot('Data', 'Precipitacao', df, 'Precipitação em mm')
Vamos agora entender o consumo de cerveja:
line_plot('Data', 'Consumo_cerveja', df, title = 'Consumo de cerveja em 2015')
corr = df.corr()
corr
# matriz de correlação de Pearson
correlation_matrix = df.corr()
mask = np.zeros_like(correlation_matrix)
mask[np.triu_indices_from(mask)] = True
sns.heatmap(
correlation_matrix,
annot = True,
mask = mask,
square = True,
cmap = 'Blues');
df.iloc[:,1:].plot(kind = 'box', figsize = (12, 5));
df.iloc[:,1:].hist(figsize = (12, 8), bins = 30);
sns.pairplot(df);
Criando o modelo¶
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
X = df.drop(columns = ['Data', 'Consumo_cerveja'])
y = df['Consumo_cerveja']
X_train, X_test, y_train, y_test = train_test_split(X, y,
test_size = 0.3,
random_state = 3)
X_train.shape, y_train.shape, X_test.shape, y_test.shape
linear_model = LinearRegression()
linear_model.fit(X, y)
Coeficientes estimados:
predict = linear_model.predict(X_test)
Score:
linear_model.score(X_test, y_test)
mean_squared_error(y_true = y_test, y_pred = predict)
Usando statsmodels¶
import statsmodels.formula.api as smf
Vamos estimar o seguinte modelo:
$$ \text{consumo_cerveja} = \beta_0 + \beta_1\text{temperatura_mean}+ \beta_2\text{temperatura_min}+\beta_3\text{temperatura_max}+\beta_4\text{final_semana}+\beta_5\text{precipitação} $$# modelo com constante
model_1 = smf.ols(
formula = """
Consumo_cerveja~Temperatura_mean+Temperatura_min+Temperatura_Max+Precipitacao+Final_semana
""",
data = df
).fit()
model_1.summary(title = 'Modelo 1')
Agora vamos estimar o seguinte modelo:
$$ \text{consumo_cerveja} =\beta_1\text{temperatura_mean}+ \beta_2\text{temperatura_min}+\beta_3\text{temperatura_max}+\beta_4\text{final_semana}$$model_2 = smf.ols(
formula = """
Consumo_cerveja~0+Temperatura_mean+Temperatura_min+Temperatura_Max+Precipitacao+Final_semana
""",
data = df
).fit()
model_2.summary(title='Modelo 2')
- Os parâmetros do modelo são significativos.
- Log-Likelihood: usamos para fazer comparação entre modelos; varia de $-\infty,+\infty$ e quanto maior o valor, melhor o modelo. O modelo 1 é melhor que o 2 sob essa regra.
- AIC e BIC são utilizados para comparar modelos, com fundamentação estatística mais rigorosa que o $R^2$ ajustado. O modelo com menor critério de informação é o melhor.
- Teste Omnibus: normalidade da distribuição dos resíduos do modelo. 0 indica normalidade perfeita. Prob(Omnibus) mensura a probabilidade de os resíduos serem normalmente distribuídos. 1 indica distribuição perfeitamente normal.
- Skew: mede o grau de simetria (assimetria) da distribuição dos resíduos. Se a distribuição for normal a assimetria é nula. Se a distribuição for assimétrica positiva (negativa) os valores são positivos (negativos).
- Curtose: mede o grau de achatamento da distribuição de probabilidade.
- Os resíduos do modelo são normalmente distribuídos?
- Teste de autocorrelação DW:
- Valor próximo de 4: há evidência de autocorrelação negativa;
- Valor próximo de 2: há evidência de ausência de autocorrelação;
- Valor próximo de 0: há evidência de autocorrelação positiva.
- Teste F: Teste de significância conjunta dos parâmetros.
- $H_0: \beta_1=\beta_2=\beta_3=\dots=\beta_k=0$
- $H_1: \beta_k \neq 0 $
- Teste t: testa de significância individual dos parâmaetros do modelo
- $H_0: \beta_j = 0$
- $H_1: \beta_j \neq 0$
- Se $\beta_j = 0$ então não há relação linear entre a variável explicativa e a variável dependente.
predicoes = pd.DataFrame(model_1.predict(), columns=['Predições 1'])
predicoes['Predições 2'] = model_2.predict()
predicoes['Consumo de cerveja (litros)']= df['Consumo_cerveja']
predicoes.plot();
Referências
- BUSSAB, W. O.; MORETTIN, P. A. Estatística Básica. 6. ed. São Paulo : Saraiva, 2010.
- SARTORIS, A. Estatística e introdução à econometria. São Paulo: Saraiva, 2003.
- Regressão Linear: https://ivanildo-batista13.medium.com/regress%C3%A3o-linear-m%C3%BAltipla-em-python-eb4b6603a3